Force distributions in 3D granular assemblies: 
Effects of packing order and inter-particle friction 
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We present a systematic investigation of the distribution of normal forces at the boundaries of static 
packings of spheres. A new method for the efficient construction of large hexagonal-close-packed 
crystals is introduced and used to study the effect of spatial ordering on the distribution of forces. 
Under uniaxial compression we find that the form for the probability distribution of normal forces 
between particles does not depend strongly on crystallinity or inter-particle friction. In all cases the 
distribution decays exponentially at large forces and shows a plateau or possibly a small peak near 
the average force but does not tend to zero at small forces. 

PACS numbers; 8I.05.Rm, 45.70.-n, 83.70.Fn, 45.70.Cc 



I. INTRODUCTION 

In ordinary solids and confined fluids, uniformly ap- 
plied loads are distributed homogeneously throughout 
the material. This is not the case for granular materi- 
als, which are collections of discrete non-cohesive macro- 
scopic particles M. For such systems, stresses are dis- 
tributed in a highly inhomogeneous manner, along net- 
works containing the largest inter-particle forces P-pl. 
These "force chains" support most of the external load, 
effectively shielding large regions of the material. The 
probability distribution, P{F), of normal forces, F, be- 
tween neighboring particles provides a quantitative way 
of analyzing these stresses B. P{F) in granular media 
exhibits characteristic differences from what would be ex- 
pected for ordinary solids or fluids. 

Recent experiments have measured P(F) at the con- 
tainer boundaries of granular media p-M]. When the 
applied loads are small enough not to cause significant 
particle deformations two key features of P{F) have 
emerged: first, P{F) exhibits an abundance of small force 
values, even far below the mean force, F; and second, for 
forces above F, P{F) decays with a characteristic, ex- 
ponential dependence. This unusual shape of P{F) has 
wide-ranging implications. Compared to a Gaussian pro- 
file, the exponential tail indicates a significantly higher 
probability of an individual contact force greatly exceed- 
ing the mean force. In applications with a fixed yield 
strength, this corresponds to higher chances for mate- 
rial failure. Perhaps even more intriguing is the non-zero 
value for P{F) as F approaches zero. Investigations con- 
cerning the structure of force chains have suggested that 
this abundance of low forces could be caused by arching 
effects |ll|"[l4|. The shape of P{F) at low forces has also 
been associated with glassy behavior p^ ]. 

As far as the precise functional form for P{F) is con- 
cerned, however, there is currently no clear consensus. 
Several fitting forms have been proposed, which differ 



most widely in their predictions for P{F) at low forces, 
ranging from P{F) ^ [|l6| to P{F) -^ oo |ll|,0[8| 
as F approaches zero. Because of the obvious discrepan- 
cies not only among some of these model predictions, but 
also between certain model predictions and the available 
experimental data, it is important to test the robust- 
ness of the results with regard to variations in the bead 
parameters, the packing construction, and the loading 
conditions. Experimentally, the issue of whether P{F) 
depends on inter-particle friction in any significant way 
and whether disordered packings behave differently from 
highly ordered, crystalline configurations has not been 
tested. Clearly, a perfect, infinite crystal composed of 
identical grains would result in a delta function for P{F) 
since all contact forces are the same in this case. How- 
ever, it is unknown how the shape of P{F) is modified 
by the small amounts of disorder present in real granular 
crystals. Also it is not known how inter-particle friction 
infiuences the distribution of forces. 

Here we address these issues with a systematic exper- 
imental investigation of the effects on P{F) of packing 
order as well as inter-particle friction. Our results con- 
firm the robust character of the probability distribution 
of forces. We find that P{F) is essentially independent 
of changes in the particle arrangement from amorphous 
to crystalline and changes in the particle surface from 
smooth to rough. In particular, we observe an exponen- 
tial decay over more than 2.5 decades for forces larger 
than the mean. Below the mean force, our data are con- 
sistent with either a flat distribution or a small peak cen- 
tered close to F. 

The paper is organized as follows. After a review of 
currently available theoretical model predictions and sim- 
ulation results, we introduce our experimental methods, 
including details of the carbon paper technique used to 
measure normal forces and a discussion of the experi- 
mental uncertainties. We then present results for P{F) 
obtained from crystalline and amorphous granular pack- 
ings using both smooth and rough beads. We also show 



how P{F) evolves with system depth. A final section dis- 
cusses these results in light of recent experimental and 
theoretical work on granular and other jammed systems. 



II. BACKGROUND 

The original g-model, proposed by Coppersmith et 
al. plpq] captures the dominant behavior of P{F) at 
large F: P{F/F) oc exp(— dF/F), where d is a posi- 
tive constant. This model was extended by Nguyen and 
Coppersmith ||l^ and by Socolar and Sexton |^,0 , and 
Claudin et al. have related it to a diffusion equation with 
an additional randomly varying convection term 122] . De- 
spite this success, these models predict power-law behav- 
ior for P{F) at forces below the mean force, implying 
P{F) — > as F ^ 0, which disagrees with experimental 
findings [0,p|. Usingcontact dynamics computer simula- 
tions, Radjai et al. |l^jlj,|2^ have measured P{F) and fit 
it with a slowly diverging power law below F and with a 
decaying exponential at large forces. This form is quali- 
tatively more consistent with experimental results H,@| . 

Recently, the effects of particle deformations on the 
structure of force chains has been considered. Simula- 
tions by Makse et al. H show evidence for a crossover in 
the shape of P{F) from pure exponential to Gaussian as 
the applied load is increased. However, experiments by 
L0voll et al. ^ with no external force do not show pure 
exponential behavior; rather, they are more consistent 
with the results of Mueth et al. |Q at higher deforma- 
tions. This is in agreement with simulations by Thorn- 
ton |2j,g^ and Nguyen and Coppersmith [^ , which ob- 
serve a slow trend toward Gaussian behavior at very large 
deformations with no evidence of a pure exponential in 
the low-deformation regime. 

Tkachenko and Witten have argued that frictional 
forces lead to elastic behavior | 27 | . Simulations per- 
formed by Eloy and Clement |2q| and by Thornton ||2|] 
predict a dependence of P{F) on inter-particle friction, 
but this has yet to be seen in experiment. 

P{F) has also been studied in simulations of super- 
cooled liquids and glasses [^. For a hquid with a 
strongly repulsive core the tail of P{F) is a decaying ex- 
ponential at all temperatures. As the temperature of the 
liquid is lowered to the glass transition, P{F) develops a 
small peak centered at F and is very reminiscent of the 
experimental data on granular systems found by Mueth 
et al. H. 



III. EXPERIMENTAL METHODS 

We have performed experiments on five different gran- 
ular systems. In order to compare our results with the 
earlier work of Mueth et al. m, we first examined an 
amorphous pack of smooth spherical glass beads. We 
then constructed both hexagonal-close-packed (HCP), 




FIG. 1. A schematic of the triangular cell used to contain 
both crystalline and amorphous packs. A force of 7600 N was 
applied to the rigidly constrained shaft on the top piston. The 
cell is 53 bead diameters on a side. A portion of a scanned 
image of the imprints left on the white paper by a crystal pack 
is shown with a close-up highlighting the crystalline order. 



and face-centered-cubic (FCC) crystalline lattices of the 
spheres and repeated the experiment on these highly or- 
dered systems. In addition, we examined both amor- 
phous and HCP arrangements of spheres with roughened 
surfaces. 

The experiments were performed using approximately 
70,000 soda lime glass spheres with diameter 3.06 ± 
0.04 mm. We used an acrylic container of equilateral tri- 
angular cross section with a side length of 165 mrn, which 
is commensurate with our HCP crystals (See Fig. |l|). The 
top and bottom boundaries of the system were thick, 
close-fitting acrylic pistons. The height in layers, h, of 
the HCP crystals ranged from 1 to 61, with the bulk 
of the data being taken at h = 45. FCC crystals were 
studied with h = 10. The amorphous arrangements had 
a filling height of 12 cm, which corresponds to approxi- 
mately the same number of particles as in the 45-layer 
crystals. 



A. Crystal Packing 

It is difficult to construct three-dimensional macro- 
scopic crystals of granular materials. Vibration meth- 
ods for crystallization can efficiently produce triangular 
close-packed layers ^^1; however these methods do not 
control the stacking order of the planes which is neces- 
sary for three-dimensional crystal structures. Also, these 
methods typically result in large numbers of defects (ap- 
proximately 10%) which can ruin long range order. We 
have developed a method to construct near-perfect crys- 
tals of glass spheres in an efficient manner. Our tech- 
nique consists of four steps. First, one layer's worth of 
beads is placed upon a perforated metal sheet which has 
small holes arranged in a triangular lattice with a spac- 
ing equal to the bead diameter. The array of holes, and 
therefore the collection of beads, is commensurate with 
our cell. Second, a vacuum is established behind the 
perforated plate; this holds the beads in place while the 
plate is manipulated. Third, the arrangement is inverted 
and lowered into the cell, which has been precisely ma- 
chined to accept a monolayer of crystallized beads. The 
vacuum is then removed to release the beads. Fourth, 
a small number of residual defects and grain boundaries 
are removed by hand before the process is repeated to 
add successive layers. By carefully monitoring the layer- 
to-layer structure, we were able to create large, nearly 
perfect HCP structures of arbitrary height with 1,431 
and 1,378 particles in alternating layers. Upon construc- 
tion, each 45-layer HCP crystal contained 70,119 beads 
with fewer than 10 beads significantly out of place, and 
typically these defects were confined to the edges of the 
container. 

HCP packings of rough spheres were created in the 
same manner. Unfortunately, the etching process re- 
duced the bead size by a few micrometers, causing a 
noticable incommensurability between the bead lattice 
and the cell. This resulted in a larger number of defects 
and less stability. 

We also constructed a FCC crystal, but boundary 
effects made the procedure considerably more difficult. 
These effects were minimized by rotating the orientation 
of the pack by thirty degrees. A seed layer was used at 
the bottom of the cell to arrange the pack in this awk- 
ward position. The resulting packing was incommensu- 
rate with our cell and was highly unstable; for this reason 
only one small {h = 10) smooth FCC crystal was created. 



B. Rough Particles 

Force chains are present only when grains are in con- 
tact. This suggests that the highly nonlinear frictional 
forces between particles may play an important role in 
determining the network of force chains. To understand 
the magnitude of these effects on P{F), we have varied 
/i, the coefficient of static friction between grains. 




FIG. 2. Scanning electron microscope images of a) a 
smooth glass bead and b) a glass bead after etching in hy- 
drofluoric acid. The white bar is 100 ^ra in length. 



A collection of glass spheres was etched in a hydrofluo- 
ric acid solution (Ashland Chemical Inc. 7:1 Buffered Ox- 
ide Etch) for approximately 2 hours. The difference be- 
tween the smooth beads and the etched beads was quite 
apparent by visual inspection. Figure ^ shows scanning 
electron microscope images of smooth and etched beads. 
We have quantified the increase in surface friction be- 
tween grains due to etching using the apparatus shown 
schematically in Figure |3[ The ratio of the coefficients of 
static friction is characterized by the changes in spring 
length given by 
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where the change in spring length. Ax, is measured at 
the point of bead movement, and l^Xsied is included to 
remove the contribution from friction of the sled with the 
base. The measured ratio was 3.0 ± 0.3 . 
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FIG. 3. Schematic of the apparatus for measuring static 
friction between two beads. Fn is the normal force applied 
to the beads, and Ax denotes the spring elongation required 
to move the bottom bead which is attached to a low-friction 
sled. 



C. Analysis Techniques for Measuring Local Forces 

The experiments were performed as follows. An ex- 
ternal force of 7600 N was applied to the top piston of 
the cell with a pneumatic press, while the bottom piston 
rested upon a fixed floor. Some force was carried by the 
walls of the cell, but we determined that this was typi- 
cally less than 10% of the applied force. We measured 
the normal forces at the upper and lower boundaries by 
placing carbon paper and white paper between the gran- 
ular material and the pistons [3^]. This allowed each 
grain at the boundaries to press into the paper, leaving 
a mark whose size and intensity were dependent upon 
the normal force on that bead P,M,p2]. A portion of a 
scanned image of the resulting imprints is shown in Fig- 
ure E, Roughly 2000 marks from both the top and bot- 
tom boundaries were identified for each data run. We 
observed no significant deviations from circular marks, 
as might be produced by tangential forces. 

After each data run, each sheet of white paper was 
carefully removed and digitized with a flat bed scanner. 
These digital images were then processed using image 
analysis software to find both the area and the inten- 
sity of each mark. The force on each bead was extracted 
from the area and intensity of the corresponding feature 
using calibration curves taken over the same force range. 
To account for variations in the applied load from run 
to run, the forces on the individual beads were divided 
by the average force for each data run. A histogram of 
the rescaled forces was then averaged over typically 10-20 
independent data runs and normalized in order to rep- 
resent the probability distribution of forces, P{f), where 
/ = F/F. ^ 

Calibration data were obtained by lowering a known 
mass connected to a single bead onto carbon paper 
covering white paper. This system, described in ref- 
erence 0], consisted of a vertical slide with impulse- 
absorbing springs; a pulley and cable system was added 
to ensure that the bead was lowered slowly and mono- 
tonically in order to further eliminate impulse forces. Be- 
cause of the potential variability in carbon paper, each 
package of paper was calibrated independently. To inter- 
polate a calibration curve from the discrete set of calibra- 
tion points we used two procedures: 1) a global fourth- 
order polynomial fit and 2) a local sliding parabolic fit 
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FIG. 4. Calibration data, a) A typical calibration set fit to 
a global fourth-order polynomial (GFP). b) Percent difference 
between the fit in the top panel and a local sliding parabolic 
fit (LSP) to the same data, given by 100 (LSP - GFP)/GFP. 
While both interpolations fit the data reasonably well, the 
differences below the mean force (F ~ 7N) shown in the 
bottom panel can result in large changes to the calculated 
probability distribution. 



with a window of eleven data points. Figure ^ shows a 
representative set of calibration data together with the 
fitted fourth-order polynomial. The relationship between 
force and intensity is monotonic but nonlinear. We found 
that the force distributions were very sensitive to the cal- 
ibrations, especially in the low-force regime. 



D. Systematic and Statistical Errors 

The uncertainties due to the calibrations are an inher- 
ent characteristic of the carbon paper method and rep- 
resent the largest contribution to systematic errors. The 
most significant problem is the proper interpolation be- 
tween calibration points. As shown in figure ^b, the two 
interpolation methods can lead to significant differences 
in the intensity-to-force conversion. It is not clear which 
of the two fitting methods gives the more accurate cali- 
brations, although we found the fourth order polynomial 
to produce smoother data. For the data shown below, we 
computed P{f) using both methods (see insets of figures 
H and o), and we interpret the difference as an estimate 
of the systematic error introduced by the calibration. 



An additional complication arises from the presence 
of forces below our threshold of detectability, which is 
approximately 0.7 N. It is difficult to distinguish between 
undetectably small forces and missing bead contacts, but 
the number of forces below threshold can be estimated 
experimentally using the two-sided tape technique de- 
scribed by Mueth et al. p|. An approximate number of 
very small forces was inserted into the force distributions 
to account for these undetectable contacts. Errors in 
the number of undetectable forces affect the normaliza- 
tion and thus effectively rescale /. Including this source 
of error, we estimate the uncertainty in the exponential 
slope to be approximately 5%. 



IV. RESULTS 

We assembled a large number of amorphous and crys- 
talline packs of both smooth and rough glass beads and 
studied them with the carbon paper technique. As de- 
tailed below, the force distributions exhibit exponential 
behavior for / > 1. This result has been found previ- 
ously in experiment, theory, and simulation ||7|, p^ , p^ , p3[ . 
Our data are consistent with a small peak in P{f) near 
the mean force or with a plateau below the mean force. 
We found that the force distributions, P{f), were well 
represented in all cases by the functional form 



P{f)^a{\-he-'i")e 



-df 
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which is a slight generalization of that proposed by 
Mueth et al. [§. 

Data was taken at both the top and bottom boundaries 
of the pack. Very little force was carried by the walls, so 
that the total force at the bottom, as measured from 
the carbon dots, was typically within 15% of the force 
at the top. As had been found previously 0, P{f) was 
equivalent between the two surfaces, and the following 
data represent an average of the results at the top and 
bottom in order to improve the statistics. 

The total forces on both boundaries were typically 
within 15% of the applied force. This result implies that 
the carbon marks were dominated by normal forces, as 
additional tangential forces would have increased the to- 
tal force. The equivalence between total measured force 
and applied force also suggests that there were not a sig- 
nificant number of duplicate marks resulting from bead 
movement as these would also tend to raise the measured 
total force. 



A. P(/): Amorphous vs. Crystalline 



The data shown represent an average over 32 indepen- 
dent packings. Each run gave approximately 1800 force 
measurements. The distribution agrees well with that re- 
ported by Mueth et al. . The solid line in the figure is 
a fit to Equation 0; the fitting parameters arc shown in 
Table |. The error bars represent the statistical error for 
multiple data sets, and the insets show P{f) for low forces 
calculated using both calibration interpolation methods. 
We also measured the correlations between pairs of nor- 
mal forces, and we found no significant correlations for 
amorphous arrangements. 

P{f) for multiple HCP crystals, h — 45, of smooth 
beads is shown in figure 0b. This represents an average 
over 26 runs, distributed over three different HCP crys- 
tals with each run giving approximately 2200 force mea- 
surements. These data show a similar form to that of 
the amorphous packs, but there appears to be a slightly 
higher peak in the smooth crystal data for forces near 
the average value, / w 1, as well as a slightly steeper 
exponential decay at large /. We have seen some evi- 
dence of aging (a gradual settling of the grains into a 
more jammed state) for a single HCP crystal which had 
been previously subjected to more than 40 runs. We have 
excluded these data from our results. 

To test the robust nature of the probability distribu- 
tion, we constructed an FCC crystal with h ~ 10. Fig- 
ure oc shows P{f) to be virtually indistinguishable from 
that obtained for the HCP crystal. Due to the difficul- 
ties of constructing a crystal incommensurate with our 
cell, only one FCC crystal was created. The data repre- 
sent an average of 10 runs with approximately 2000 force 
measurements each. 



B. P{f): Smooth vs. Rough 

Amorphous and HCP crystalline packs created with 
the roughened beads were also studied. The rough 
amorphous data were taken from 20 independent pack- 
ings, each representing approximately 1800 force mea- 
surements. The rough HCP crystal data were taken 
from 10 data runs on one crystal with approximately 
2200 force measurements each. The resulting force dis- 
tributions for the rough beads (Fig. ph were not signif- 
icantly different from the distributions for the smooth 
beads (Fig. ||). The fitting coefficients are given in Ta- 
ble |. 

All systems, consisting of both smooth and rough 
beads in both amorphous and crystal packs, showed sim- 
ilar force distributions. The P{J) curves for four systems 
are shown together in Figure \ 



Amorphous packs were recreated before each run, 
which ensured that each run was independent of the oth- 
ers. Figure |p, shows P{f) for smooth amorphous packs 
which have been calibrated using the polynomial method. 



C. P{f): Height Dependence 

The aforementioned experimental results were ob- 
tained for HCP crystals and amorphous packs which were 
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FIG. 5. P{f) for a) amorphous packs of smooth glass 
spheres, b) HCP crystalhne packs, h = 45, of smooth glass 
spheres, and c) an FCC crystalline pack, h — 10, of smooth 
glass spheres. The error bars represent the statistical uncer- 
tainty, and the solid line represents a fit to Eq. bl Insets show 
the difference between the fourth-order polynomial (o) and 
sliding parabolic (•) interpolation methods. 
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FIG. 6. P{f) for a) amorphous packs of rough glass spheres 
and b) HCP crystalline packs, h — 45, of rough glass spheres. 
The error bars represent the statistical uncertainty, and the 
solid line represents a fit to Eq. y. Insets show the dif- 
ference between the fourth-order polynomial (o) and sliding 
parabolic (•) interpolation methods. 



Fitting Parameters 



Smooth Amorphous 


1.5 


0.59 


3.1 


1.21 


Smooth HCP Crystal 


2.7 


0.80 


2.0 


1.50 


Smooth FCC Crystal 


2.8 


0.80 


1.5 


1.48 


Rough Amorphous 


2.4 


0.75 


1.4 


1.41 


Rough HCP Crystal 


3.9 


0.85 


0.8 


1.60 



TABLE I. Values of fit parameters for the form 
P{f) = a [1 - fe exp(-c/2)] exp(-d/) 
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FIG. 7. Force distributions for smooth amorphous (o), 
smooth HOP (•), rough amorphous (a), and rough HOP (a) 
granular packs. 



many bead diameters in height. We are interested in un- 
derstanding how the P{f) distributions are built up as 
the height is increased. We studied this question by vary- 
ing the number of layers in HCP crystals of smooth glass 
beads. Figure || shows P(/) for HCP crystals of different 
heights ranging from h — 1 to h = 61. For comparison, 
the solid line in each panel represents the fit to Equations 
for the /i = 45 smooth HCP crystals. For a single layer, 
h = I, P{f) is determined primarily by the polydisper- 
sity of the beads themselves. As the height is increased, 
P(/) evolves toward a height-independent form. We find 
that P{f) is robust for pack heights greater than h « 15. 




0123456701234567 
/ / 

FIG. 8. Force distributions, P{f), for smooth HCP crystals 
of various heights, h. The solid line is a fit to Eq. ^ for the 
/i = 45 smooth HCP crystal. 



D. Particle Deformation 



Recent simulations on both 2-D and 3-D bead packs 
suggest that the form of P{f) depends on particle de- 
formations ||9|,p4|-p6|| . We calculate particle deformations 
using a Hertzian deformation law |33|: 



1.53 mm, and F « 7.0 N, resulting in particle deforma- 
tions of 0.2%. We varied the force applied to the h — A5 
HCP crystals over a range of F from 3.5 N to 8.0 N. This 
corresponded to deformations of 0.12% to 0.2%. Over 
this limited range, we found that P{f) remained un- 
changed. 



5R 



2D^F 



where F is the mean force per bead, and R is the bead 
radius. D is determined by material properties: 



D 



3(1 -g^) 
2E 



where a is Poisson's ratio and E is Young's modu- 
lus. Soda lime glass spheres have a = 0.24 and E — 
69GPa B4|. For the data we have presented, R = 



V. DISCUSSION AND CONCLUSIONS 

Our data for P{f) can be well fit by the form given 
in Equation 0. For forces larger than the mean force, 
we found P{f) to be a decaying exponential. For forces 
smaller than the mean force, our data are consistent with 
either a small peak or a plateau. There was very little 
change in the force distribution when the bead pack was 
varied from amorphous to crystalline, although there is 
some evidence that the exponential decay for large forces 
is slightly steeper in the crystalline configurations. There 



may also be support for an aging effect, although this is 
still inconclusive. We found P{f) to be unaffected when 
the inter-particle coefficient of friction was changed by a 
factor of three. 

Packing history may also affect the results for an in- 
dividual data run and could therefore alter the form of 
P(/). A small number of particular runs were found to 
show significantly different results, such as pure exponen- 
tial behavior. Only after a considerable number of indi- 
vidual experiments were performed and averaged could 
a robust probability distribution be presented. We find 
these variations in the individual runs to be quite surpris- 
ing. They may indicate that the system is occasionally 
caught in a state not representative of the average con- 
figuration. 

We do not see evidence for a change in P{f) due to the 
deformation of the particles over a range from 0.12% to 
0.2%. Simulations and experiments performed by Makse 
et al. suggest that P{f) is a pure exponential for small 
deformations, with a transition to a Gaussian centered 
at / = 1 when particle deformation is increased [g|. We 
estimate the deformations to vary from 0.09% to 0.4% 
in the experiments of Makse et al., while their simula- 
tions have a much wider range. The fact that our data 
do not show a transition toward either pure exponen- 
tial or Gaussian behavior seems to be in contradiction 
to their experimental results. Our data are not neces- 
sarily inconsistent with their simulations, as our results 
may be viewed as being in a transition regime between 
small and large deformations. We note, however, that 
simulations by Thornton p4,p5[ as well as Nguyen and 
Coppersmith p6| find a slow transition toward a Gaus- 
sian profile at high deformations (well beyond our range) , 
but do not find a pure exponential in any deformation 
regime. At low deformations, these simulations show a 
form of P(/) which is qualitatively similar to our results. 
In addition, this form is also consistent with experiments 
by L0voll et al. S on a granular pack with no external 
forcing and thus very low deformations. 

We are able to make qualitative comparisons with the 
predictions of some theoretical and simulation work for 
P(/) in the / < 1 regime. We find that P{f) does not 
approach zero as / -^ 0, in contrast with the predic- 
tions of the g-model flH]. Based on simulations, Radjai 
et al. report a slow power-law divergence as / approaches 
zero [0,18U23]; while this form is similar to our exper- 
imental results (since the power-law exponent is very 
small), it does not account for the possibility of a peaked 
distribution near the mean. 

Simulations performed by O'Hern et al. ||lj] on 
quenched molecular liquids in two dimensions show the 
most apparent similarities to our results. For potentials 
with finite repulsive terms, such as a Lennard- Jones po- 
tential with a range cut off at its minimum value, the 
P{f) distribution is indistinguishable from our findings. 
Also, there may exist a correlation between granular crys- 
tals which have aged (as mentioned above) and quenched 
molecular liquids at temperatures well below the glass- 



transition temperature. The increase of the exponen- 
tial slope as well as a slight peak near the mean force 
in crystalline configurations are similar to those found 
in jammed fluids where large forces cannot relax when 
quenched below the glass-transition temperature. 
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